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APPLICATION OF THE METHOD OF INTEGRAL RELATIONS TO 


SUPERSONIC NONEQUILIBRIUM FLOW 

PAST WEDGES AND CONES 

By Jerry C. South, Jr. 
Langley Research Center 


SUMMARY 


The method of integral relations is used to calculate supersonic nonequi- 
librium flow past wedges and circular cones at zero incidence. Vibrational 
relaxation in a pure diatomic gas is considered so that the results can be com- 
pared with existing calculations using the method of characteristics. The 
approach taken is similar to that followed previously in ideal-gas problems, 
except that the isentropic law cannot be used and the vibrational rate equation 
is included. The governing equations are converted to approximate systems of 
ordinary differential equations which are solved as initial-value problems on a 
high-speed digital computer. Alternate procedures are discussed which point 
out certain features of general importance. 

Numerical results are presented to illustrate the convergence and accuracy 
of the method in predicting the distributions of flow variables during the 
approach to equilibrium. The calculations are performed through the third 
approximation for the wedge and through the second approximation for the cone. 


INTRODUCTION 


The method of integral relations (ref. l) is a numerical analysis technique 
for solving the nonlinear equations of gas dynamics. Specifically designed for 
high-speed computing, the method has been shown to be useful and versatile in a 
variety of problems. Most applications have been in mixed flows, such as the 
sonic flow past ellipses and ellipsoids (ref. 2) and the supersonic blunt-body 
problem (refs. 3 to 9)- Other important applications have been made to conical 
flows without axial symmetry (ref. 10 ) and to the viscous boundary layer 
(ref. 11). 

In the aforementioned works the gas was considered to be ideal; thus, an 
obvious step would seem to be the extension of this method to include real-gas 
effects in hypersonic flow problems. There do not appear to be any obstacles 
to such an extension in the case of equilibrium thermodynamics. The entropy is 
still conserved along the stream l ines and, therefore, the stream function can 



be used advantageously (refs. 3 and 6). Usually the variations of equilibrium 
flow properties across the shock layer are no' more nonlinear than are those in 
ideal-gas cases. In nonequilibrium flows , however, there are some difficulties: 
the isentropic law is not valid and cannot be used; chemical and vibrational 
rate equations must be included; and shock-layer profiles can be highly non- 
linear. This last difficulty is particularly important to a basic feature of 
the method: the assumption of interpolation polynomials for certain groupings 

of the flow variables. The nonlinear profiles may be poorly approximated by 
low-degree polynomials and might considerably retard the convergence of the 
method. On the other hand, the integral relations effect a "smoothing” of pro- 
file irregularities (ref. 11); thus it seems possible that some significant 
results might be obtained in low-degree approximations. 

Two applications of the integral relations to nonequilibrium flow calcula- 
tions have recently appeared. Shih et al. (ref. 12) used the first approxima- 
tion to calculate hypersonic flow of a five -component (N2,02,N,0, and NO) dis- 
sociating gas past a sphere. Some of the different procedures necessary for 
nonequilibrium calculations, as compared with previous ideal-gas work (refs. 3 
to 9) , were emphasized, particularly in regard to the starting conditions at the 
axis of symmetry. 

South (ref. 13) applied the integral relations to supersonic, nonequilibrium 
flow past wedges and cones, where the molecular vibrations are relaxing. As in 
reference 12, only the first approximation was considered, but the results for 
shock-wave shape and surface pressure distribution compared favorably with cal- 
culations based on the method of characteristics (supplied by R. Sedney and 
N. Gerber of the Ballistic Research Laboratories, Aberdeen Proving Ground). In 
reference 13 the vibrational relaxation time was assumed to be constant through- 
out the wedge or cone shock layer to eliminate dependence on experimental results 
for the relaxation time. The computational program was later modified to account 
for a pressure- and temperature-dependent relaxation time, and the first approxi- 
mation results were found to be adversely affected. Similar unpublished calcu- 
lations made at the Langley Research Center using the Li ghthi 11 -Freeman dissoci- 
ating gas model showed that the first approximation was inadequate. It appears 
that in many problems approximations higher than the first may be needed. 

In the present paper, the work of reference 13 is extended to higher approx- 
imations in a straightforward manner. The prime objective is to observe the 
resulting gains in overall accuracy, with emphasis on certain important features 
of the method of integral relations that have not been evident previously. 


SYMBOLS 


Cp frozen-flow specific heat at constant pressure 

E vibrational energy 

E e q equilibrium vibrational energy 
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J 

M 

N 

P 

Q 

R 

r 

T 

u 7 v 

V 

x,y 

Ax 

P 

7 

6 


nonhomogeneous functions in equations (l) to ( 4) , used, in equations (l4) 

to ( 19 ) 

functions differentiated with respect to y in equations (l) to (4), 
used in equations (l4) to ( 19 ) 

0 or 1 for a wedge or a cone , respectively 
frozen-flow Mach number. 


fr p/p 

order of approximation (number of strips) 
pressure 

functions differentiated with respect to x in equations (l) to (4), 
used in equations (l4) to (19) 

gas constant 

radial coordinate normal to cone axis 
temperature 

velocity component in x- and y-direction, respectively 


total velocity , \j u^ + 

coordinate along and normal to body surface, respectively 
integration step size 
shock-wave angle 

ratio of frozen-flow specific heats 
shock-layer thickness in y-direction 


E ea - E 

vibrational driving force, — 


© v characteristic vibrational temperature 

0 wedge or cone half -angle 

A shock-layer included angle, p - 0 

p density 
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T 


vibrational relaxation time 


\jr stream function 

Subscripts : 

00 free -stream quantity 

1 1,2, 3,4 (refers to continuity, x-momentum, y-momentum, and rate equa- 

tions, respectively) 

0 surface quantity; see definition of subscript k 

k 0,1,. . . ,N-1,6 (denotes a particular strip boundary as labeled in 

fig. 2) 

5 shock-wave quantity; see definition of subscript k 

Primed quantities are dimensional; unprimed quantities are dimensionless, 
as shown after equation (6). 

Barred quantities are corrected surface flow variables, 

PROBLEM DESCRIPTION 


The physical problem to be studied is the steady supersonic flow of a pure 
diatomic gas past symmetric wedges and right circular cones. The flow is 
inviscid and the only dissipative mechanism is the relaxation process of the 
molecular vibrations. This is then the same problem studied in references lh 
and 15 by using the method of characteristics. There are several reasons why 
this model is well suited to the objectives of the present study. The governing 
equations are relatively simple in form, yet the physical problem is highly non- 
linear and exhibits many of the important features of more complex none qui librium 
flows. The added difficulties which arise in the blunt-body problem (refs. 3 to 
9 and 12), such as unknown initial conditions and multiple singular points, are 
not encountered here. This is a great advantage in numerical calculations, 
particularly in higher approximations. Finally, it seems worthwhile to point 
out that radiation heating problems (refs. 16 and 17) have caused a renewed 
interest in pointed configurations, such as the wedge or the cone, for the 
design of hypersonic vehicles. 


Basic Equations 

The flow geometry and coordinate system are illustrated in figure 1 for a 
cone. A body-oriented coordinate system, with x and y the coordinates along 
and normal to the body surface, is used so that the X-axis is inclined to the 
stream direction by the angle 8 and coincides with the wedge or cone surface. 
The wedge or cone tip lies at the point (0,0). The basic equations are as 
follows (ref. 13): 
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Continuity: 


sH + i( pvrJ ) = 0 


(i) 


x -momentum: 


bx 


(p + pu 2 )r*|j + J^puvr ^ j 


jp sin 0=0 


( 2 ) 


y -momentum : 


J^puvr^ + ~ ^p + pv 2 ^rJ - j 


jp cos 0=0 


(3) 


Vibrational rate: 


^-^puEr^ + •^■^pvEr^ - per^ = 0 


(4) 


Energy: 


T + E + 7 1 Moo 2 (u 2 + V 2 ) = 1 + 7 g 1 Me,, 2 


(5) 


State : 


7M C0 2 p = pT 


(6) 


where j = 0 or 1 for a wedge or cone, 
respectively, and r = x sin 0 + y cos 0. 
The variables and coordinates have been 
nondimensionalized as follows: 


u' ,v T 


U , V = Z-L 


P = 


iL 


p» ,v « 


T = 


... T' 


E = 


E' 


I _ t 

c n 


x,y 


p = -4 


= x, ;y' 


The primes denote dimensional 
quantities. The length scale L 1 is 
taken to be V 0o t T , (0,0) > where r f (0,0) 

is the vibrational relaxation time on 
the surface at the wedge or cone tip. 

It is assumed, as in references 14 
and 15 y that 


P't 1 oc exp(C ! /T l ) 


! /ml 


(7) 


Shock wave-^ 


Streamline 



Figure 1. - Flow geometry and coordinate 
system for a cone. 
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where the constant C 1 is determined from experimental data. The vibrational 
"driving force" e appearing in equation (4) is (as in refs. 14- and 15 ) 

Epn — E 

6 = -Si (8 

T 

where t = t t /t ? (0,0) and E e q is the vibrational energy that the flow would 
have locally if it were in equilibrium at temperature T and is taken to be 



(9) 


The factor 2/7 is the ratio R'y/cp 1 


for an ideal diatomic gas (7 = 7/5)* 


Boundary Conditions 

It is assumed that the molecular vibrations are frozen across the shock 
wave and that E OT = 0; then, at the shock curve, or y = 6 (x), 

E(x,5) = Eg(x) = 0 (10) 

The ideal-gas shock-wave relations apply and the other variables can be deter- 
mined as algebraic functions of Moo, 0, and the shock-wave angle p(x). The 
necessary relations are listed in appendix A. 

The location of the shock curve is determined by 


— = tan A (11) 

dx 


where A = 0 - 0, and the shock wave is attached at the body tip so that 


5(0) = 0 


(12) 


The wedge or cone surface is a streamline; thus 


v(x, 0 ) = v 0 (x) = 0 


(13) 
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APPROXIMATE SYSTEMS 


The method of integral relations provides a logical procedure for con- 


verting equations 


Shock wave-.. 



Figure 2.- Nth-strip arrangement. 


approximate systems of ordinary differential 
equations, which are readily adaptable to 
electronic digital computing. The conversion 
procedure for the present problem is the same 
as that used in previous ideal-gas, mixed-flow 
studies (refs. 1 and 3)* The first approxi- 
mation was derived in detail in reference 13 , 
and higher approximations used herein are 
straightforward extensions of that work. 

Equations (l) to (k) are in the diver- 
gence form - that is. 




dx 




Fj_ = 0 


(i = 1,2, 5, 4) 


( 1*0 


where the appropriate functions , G^, and 

Fj_ can be seen by inspection of equations (l) 
to (4). In the Nth approximation the region 
of flow between the shock wave and body sur- 
face is divided into N equal strips, 
arranged as in figure 2. Each of equa- 
tions (1*0 is integrated from the surface 
(y = 0) to the boundary of each strip. The 
result is (for the present problem) 4N inte- 
gral relations, as follows: 


d r^k r y k 

+ - / F, dy = 0 


( 15 ) 


(i = 1,2,3, 4) 

(k = 1,2,. . . ,N-1,5) 

To evaluate the unknown integrals. Nth-degree interpolation polynomials in y 
are assumed for the functions Qj_ and Fj_j for example. 


N 

Qi = 51 < 11 , 11 ^ 

n=0 


(16) 


where the coefficients q^ n depend linearly on the x-dependent functions Qi^k* 
An approximate system of 4 n ordinary differential equations is finally obtained, 
where x is the independent variable. 
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Corresponding to the general form indicated by equations (l4), the 
three approximations (W = 1,2,3) can be written as follows: 

N = 1: (i = 1,2,3, 4) 

8 + 8 IbT^ + ( 4i '° ' 4i > B )S ' 2 ( Gi <° ' Gl ’ 8 ) ‘ 6 ( Fl >° + Fi ' 5 

N = 2: (i = 1,2,3, 4) 


N = 3: 


- **(<>! )0 - “iA + <> 1 ,b) - 8 ( F 1,0 - F i,s) * 0 

28 ^ +8 ^ + ^-<m)I 

- G±,0 “ ^ G i,l + 5Gi^ 5 - 5^2F i# i + Fi,g) = 0 

(i = 1,2,3, 4) 

28 ^ + 26 ^ + ( 2Q i,o - 9%,1 1 18%,2 - u %,s)i 

- 13 G i,0 + 270!, 1 - 27Gi,2 + 13Gi , 8 - 28^,0 + F ± , 5 ) = 0 

58 ^ ' 8 ^ + 6 («1,1 - 2 «i,2 + Mi 

- 2G i ^ 0 - 9G± f 2. + lSd± f 2 ~ 7G^ 5 - 8^3F ifl - = 0 

6S ^2 + 28 + 5(%, 1 + ^, 2 - 5% (8 )f 

+ G i,0 - 9Gi,l " 9Gi^ 2 + 17 G i , 5 “ 26 (3Fi,2 + F i,s) = 0 


first 


= 0 

(17) 


(l8a) 


(l8b) 


( 19 a) 


( 19 b) 


( 19 c) 
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The second subscript refers to a particular strip boundary, as indicated in 
figure 2. Since the shock-wave functions Q^,g are functions of p (and rg 
when j = 1 ), equations ( 17 ) to ( 19 ) were manipulated so that derivatives of 
Q-j^g and only one other function, Qi,k (k - 0,1,* . . ,N-l), appear in each 

equation. This form is algebraically convenient because dp/dx can be obtained 
easily from the single equation containing c ^3,o/ dx and dQ^g^dx (e.g, , 

eq. ( 17 ), (l8a) , or ( 19 a)) by noting that Q^o = Pq u o v o = 0 for all values 

of x. The algebra required to obtain the remaining derivatives of u^, v^, 

etc., is considerably simplified. Further details are given in appendix B, 
where the final computational equations are derived for the second approximation, 
N = 2. 


It should be noted that in the approximate systems the functions Q-j: ^ 
contain the factor r^ . If the derivatives of r k J are expanded and the equa- 
tions are then divided by r qJ , the resulting differential equations for the 

cone are very similar to those for the wedge, with some added terms (for example, 
see ref, 13 for N = 1 and appendix B for N = 2). 


Tip Solution 

The condition that the shock wave is attached, or 5(0) = 0, causes the 
coefficients of the derivatives to vanish at the tip of the wedge or cone. If 
a regular solution exists, it is necessary that the remaining terms also vanish 
at x = 0. This condition yields 4-N algebraic equations in the 5N unknowns: 
^ 0 * • • •? U N-1* v i* • • • j V N-1* 3?0>* * • > PN-1* Po 9 ' " Pxj-i> an d 

Eq,. • • > %-l v 5^ e_t c., are considered to be functions of p through 

the shock-wave relations, with Mb and 0 specified). Eliminating T^ between 
equations ( 5 ) and (6) on each strip boundary gives N additional equations; this 
step completes the initial solution. 

The N algebraic equations which correspond to the rate equation (i = k) 
together with equation (10) yield = 0; that is, the flow is frozen through- 
out the shock layer at the wedge or cone tip. The remaining equations yield the 
exact solution for frozen flow over a wedge when j = 0; that is, u^ = ug, 

v k = = °> Pk ” Pg> e ' tc - When j = 1, an approximate solution for frozen 

flow past a cone is obtained, and this solution is especially interesting. 

In the classical problem of supersonic frozen or equilibrium flow past an 
axi symmetric cone, similarity considerations show that the shock wave is straight 
and that flow properties are constant along rays emanating from the tip. Appli- 
cation of the same considerations to the approximate systems for the continuity 
and momentum equations is straightforward, since all strip boundaries are 
straight rays when the shock is straight. Discarding the x-derivatives of the 
flow variables along the strip boundaries yields the same solution for all 
values of x as the limit solution (x 0) mentioned previously. This 
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particular application of the integral relations was also pointed out by Chushkin 
and Shchennikov (ref. 10). It is worth noting that where the exact classical 
solutions result from a boundary-value problem for nonlinear ordinary differen- 
tial equations, the integral relations yield nonlinear algebraic equations which 
are solved by trial and error. Some results of these solutions are shown in 
figure 3, wherein the first (N = 1) and second (N = 2) approximations for frozen 
flow over a cone are compared with the exact solution (ref. 18) . In reference 19 
it was shown that the first approximation agrees well with the exact solution 
(ref. 20 ) for full equilibrium airflow over a cone. 



Figure 3.- Shook-wave angle as a function of cone half -angle for frozen flow. 

For the present problem in which the entire shock layer is supersonic, 1 
the initial values of the flow variables are completely determined by the approx- 
imate systems and the shock-wave relations. Since the original partial differ- 
ential equations are of hyperbolic type (ref. 22), the integral relations prop- 
erly form an initial -value problem. (In the blunt-body problem, with subsonic 
flow in the nose region and supersonic flow on the afterbody, the partial differ- 
ential equations are of mixed (elliptic and hyperbolic) type; and the integral 
relations thus form a boundary-value problem (refs. 3 to 8).) 

^At a given value of Moo > 1, there exists a limited region of wedge or 
cone angles which are smaller than the detachment angle, yet subsonic flow 
occurs in the shock layer (ref. 21, p. 683 ). This situation is excluded in the 
present formulation. 
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The initial derivatives at x = 0 can he evaluated by differentiating the 
approximate systems once more and taking the limit (x ->0). The initial deriv- 
atives so obtained for the wedge are exact in all approximations. This result 
is not surprising, since the frozen-flow variables at the wedge tip are constant 
in the interval 0 S Y ^ S (and 5 0); thus, the interpolation polynomials 

assumed for pu, puv, etc., are exact at x = 0- The wedge initial derivatives 
on any strip boundary are derived in appendix C. 

The present method also yields linear algebraic equations for the cone tip 
derivatives. The approximate expressions are more lengthy than those for the 
wedge and are not given herein; they are given in reference 13 for N = 1. The 
exact solution for the cone tip derivatives can be obtained in a manner similar 
to that for the wedge (see appendix C) by solving a boundary-value problem for 
linear ordinary differential equations. The integration must be performed 
numerically, however, since the variable coefficients of the differential equa- 
tions depend on the nonlinear conical-flow solution. The present approximate 
calculations for the cone tip derivatives for N = 2 were found to be in excel- 
lent agreement with exact numerical solutions recently obtained by R. Sedney and 
N. Gerber of the Ballistic Research Laboratories, Aberdeen Proving Ground. 


Nonconvergence of the Rate Equation 

At large distances from the wedge or cone tip, it can be intuitively rea- 
soned that throughout most of the shock layer the flow variables should approach 
their equilibrium values - that is, those values which arise from assuming zero 
relaxation time behind the shock wave. For a finite relaxation time, there are 
three zones which appear in the flow: (l) a relaxation zone just behind the 

shock wave, (2) an equilibrium zone between the shock wave and surface, and 
(3) an entropy layer adjacent to the surface (ref. Ik) . Far from the tip, the 
equilibrium zone (2) will occupy most of the shock layer so that by comparison 
the relative thicknesses of zones (l) and (3) appear to be negligibly thin. 

Zones (l) and (3) do not actually disappear, however, and this feature causes 
difficulty in any calculation scheme. For example, the flow variables at the 
surface, excepting the pressure, do not reach the same equilibrium values as in 
zone (2), and the flow variables at the shock are always frozen even at x = «>. 
The profiles of the flow variables across the shock layer thus approach jump 
discontinuities at both the shock wave and surface as x «>. 

In the present method, this difficulty is evidenced in the approximate 
systems for the rate equation (i = k). The vibrational driving forces 

along the strip boundaries appear in the terms of the approximate sys- 

tems. Each of those driving forces should tend to zero as x increases, with 
the exception of the shock-wave driving force €5. The boundary condition at 
the shock wave (eq. (lO)) gives for all values of x 

e 5 = E . e( l? . $ > 0 (20) 

t 5 
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A closer examination of the approximate systems reveals that the coefficients 
of the physical-variable derivatives contain the factor Sr^J which is of the 
order xl + J for large values of x. In the approximate rate equations, the 

terms are a ^ so th e order x-^- + j. Then equating the coef- 

ficients of x 1+ J (that is, comparing the dominant terms for large values of x) 
leads to the following conclusion: the physical derivatives along the strip 

boundaries can decay to zero as x increases only if (k ^ 8) approach 

unrealistic nonzero values to counter the always-positive contribution of eg. 
This effect is the most severe in the surface rate equation (e.g., eq. (17), 
(l8a), or (19a)) since, as x -> 00, 


€ 0 


± 



p 0 


€ S 


or 


E 0 E eq,0 


+ 



p s T o 

p ' 0 t 5 


E eq,5 


( 21 ) 


( 22 ) 


where the plus sign corresponds to N = 1 and 3 and the minus sign corresponds 
to N = 2. Equation (22) predicts alternating "overshoots" (+) and "undershoots" 
(-) of about equal magnitude in the asymptotic value for Eq. The problem is 
less severe at the interior strip boundaries, since the F4 ^ terms (k / 0,5) 
receive more weight than the Fi^g terms. For example, in equations (19b) 
and (19c) Fij. 1 and Fl^g are weighted three times as much as Flj. 5. 


NUMERICAL RESULTS 


Three different numerical cases are presented to illustrate the accuracy 
and convergence of the method. The conditions for each case herein correspond 
exactly to the conditions for the calculated examples reported in references 14 
and 15. Thus, the results obtained by using the method of integral relations 
can be directly compared with those obtained by the method of characteristics. 
For all cases the free-stream temperature is 300° K and the gas is pure nitrogen 
(By = 11.12). The constant C = C'/Too* used in the expression for the vibra- 
tional relaxation time (eq. (7)) was obtained from the authors of references 14 
and 15. The numerical cases are as follows: 

Case I - Wedge: M «, = 6, 0 = 40.02°, C = 0.4655 x 10 ^ 

Case II - Cone: M«, = 12, 0 = 46.39°, C = 1.0137 X 10^ 

Case III - Cone: M^, = 10, 0 = 53.82°, C = 1.0137 X 10^ 
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Hyperbolic Stability Criterion 


In reference 13 only the first approximation was considered, and no numer- 
ical stability problems were encountered in cases I and II. In case III the 
initial Mach number (based on the frozen speed of sound) on the cone surface 
M o vas about 1 . 08 , and all attempts at integration away from the tip failed. 

It was believed that the instability was caused by the fact that Mq = 1 is a 
natural singular point of the system. The higher approximations were later 
discovered to be unstable not only in case III but also in cases I and II 
(Mq « 1.4 and 1.7, respectively). The problem was finally resolved by applying 

a hyperbolic stability criterion, as follows: Consider the two-strip (N = 2) 

calculation for flow past a wedge, schematically illustrated in figure 4. Three 
data points along the surface normal are 
computed at each integration step: at 

the shock, at the middle line, and on 
the surface. The left-running (C + ) and 
the right -running (C“) Mach line charac- 
teristics emanating from the three 
points intersect downstream approxi- 
mately as shown at some point (x + Ax). 

For a stable calculation, the local 
step size must not be greater than that 
given by the approximate characteristic 
mesh . For the purpose of obtaining an 
approximate criterion, it is assumed 
that all the streamlines are parallel to 
the surface and the local Mach number is 
equal to Mq. The geometry then gives 

for N = 2 

1/2 


Ax < |(mq 2 - l) ( 23 ) 


By the same reasoning it can be shown 
that for any value of N 


Ax 


< H 

- 2N\ 


Mo 


- r 


(24) 



Since the criterion given by equa- 
tion (24) is approximate, it is usually necessary to use some fraction of that 
criterion to insure success. The criterion cannot he applied at the initial 
step since S(o) = 0; therefore, a single linear step is taken at x = 0, and 
the criterion is applied thereafter. 


In figure 5 the calculated surface pressure distribution near the wedge tip 
is shown for case I and N = 2. The solid curve is an unstable calculation using 
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.588 



x 


Figure 5-- Effect of hyperbolic stability criterion. Case I. 

a first-order Euler integration procedure with a fixed step size Ax of 0.001. 
For the initial conditions of this case, that step size satisfies the N = 2 
stability criterion when x > 0.013 . Note that the solution stabilized at that 
point. The dashed curve is also a first-order Euler integration with an initial 
step size Ax of 0.001, but the N = 2 stability criterion was applied at 
every successive step, and the results are stable. 

The stability criterion was included in the computational program used in 
reference 13, and case III was successfully integrated. In reference 13 cases I 
and II were integrated (for N = 1) without the stability criterion because the 
automatic step-sizing built into the Runge-Kutta integration scheme was able to 
satisfy the N = 1 criterion within the first few steps. 


Pressure Distribution and Shock-Wave Shape 

Calculations for the wedge (case I) were performed for N = 1, 2, and 3« 
For the cone, the algebraic solutions which give the initial values and deriva- 
tives are far more tedious in higher approximations than for the wedge; there- 
fore, cone calculations for N = 3 were not attempted. 


Ik 





cases most of the approach to equilibrium 
is adequately described by N = 2, and 
the accuracy of N = 2 is considerably 
improved compared with that of N = 1. 

The nonconvergence of the surface vibra- 
tional energy, explained previously, has 
a detrimental effect on the asymptotic 
values of all the variables, but the 
effect decreases as N increases. 

Figure 9 illustrates the surface 
vibrati onal-energy di stribut ion Eq ( x ) 
for case I. The erroneous overshoots and 
undershoots mentioned earlier are clearly 
evident. 


Improvements for Certain 
Surface Variables 

In reference 13 it was pointed out 
that the approximate differential equa- 
tions are not equivalent to the correct 

x-momentum and streamline rate equations at the surface - that is, 

du Q = 1 d Pp 

dx P 0 u o ^ 

and 

ZSq = fo 

dx Uq 

The discrepancy occurs in all approximations largely because the integral rela- 
tions must account for v(Su /dy) and v(dE /Sy) throughout the shock layer. 

As was shown in reference 13, equations (25) and (26) can be used together with 
equations (5) and (6) to give improved distributions u q(x), Eq(x), Tq(x), 

and Pq(x), consistent with the pressures Pq(x) obtained from the approximate 

equations. The corrected surface vibrational-energy distribution Eq(x) is 

shown in figure and the improvement is excellent. A note of caution is nec- 
essary, however: the corrected surface variables uq, Eq, etc., are obtained 

in addition to the corresponding original variables uq, Eq, etc. It is 
tempting to replace the two approximate differential equations for uq and Eq 
by the exact equations (25) and (26). Such a "hybrid" procedure was used in the 
nonequi librium blunt-body study of reference 12, and it was tried in the early 
stages of the present work. In the present application, the hybrid procedure 
always produced an unstable system. Unbounded oscillations of the derivatives 


(25) 


(26) 
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Figure 9.- Surface vibrational-energy distribution for a wedge. Case I. 


occurred near the tip, regardless of the hyperbolic stability criterion used 
It was found that if only the approximate differential equation for E q was 
discarded and replaced by equation (26), the solutions were stable, but the 
overall results were poor compared with those obtained by using the original 
approximate equations. 


I 


x , t Coordinates 

It is known that the choice of dependent and independent variables is impor- 
tant in low approximations of the method of integral relations. If, for example, 
spherical polar coordinates are used for the cone solution (frozen or equilibrium 
flow), the results are less accurate for N = 1 than the corresponding results 
for the body-oriented x,y coordinates (see fig. 2 of ref. 13 ). 

Since streamlines have a major role in nonequilibrium flows, the x,y 
coordinates were transformed to x,\|f coordinates, where \|r is the stream func- 
tion. The integral relations were redeveloped for the wedge (j = 0) and N = 2. 
No particular advantage was evident since the strip boundaries are not stream- 
lines, and there was no difference in the numerical results for case I obtained 
by use of either the x,\|r or the x,y coordinates. 


Shock-Layer Profiles 

Even in higher approximations, the 
details of the nonequilibrium flow 
between the shock wave and the body sur- 
face are not described accurately. At 
any downstream station x there are 
only N + 1 points along a normal to 
the surface from which flow-variable 
profiles can be obtained. Likewise, the 
poor asymptotic behavior, discussed pre- 
viously, causes considerable error in 
the most interesting variable, T. In 
references 12 and 13 streamlines were 
used to obtain better details for N = 1. 

Figure 10 illustrates the profiles 
of pressure and temperature obtained in 
reference 15 for case II at a station 
x = 9*0. The three points given by the 
integral method for N = 2 are also 
shown in the figure at a station x = 8.0 
(the present calculations were not car- 
ried farther). Although the pressures 
agree reasonably well, the temperatures 
do not. The corrected surface tempera- 
ture (ref. 13) Tq is shown to agree 
well with the characteristics result. 

It appears that, if accurate 
details of the temperature and density 
profiles are required, additional cal- 
culations must be carried out along 
streamlines by using a procedure like 
that described in reference 13 . 



Figure 10.- Pressure and temperature pro- 
files for a cone. Case II. 
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CONCLUDING REMARKS 


The method of integral relations has been used to calculate supersonic, 
vibrationally relaxing flow past wedges and cones. Comparison of these calcula- 
tions with characteristics calculations showed that the second, or two-strip, 
approximation was accurate for most of the approach to equilibrium. The asymp- 
totic value of the surface vibrational energy did not converge to the correct 
value in higher approximations. The other flow variables were rather insensitive 
to this discrepancy, however, and the erroneous surface energy distribution was 
corrected. 

For the present problem, in which the entire shock layer is supersonic, 
the original partial differential equations are of hyperbolic type. The approx- 
imate systems of ordinary differential equations correctly posed an initial- 
value problem. Control of the integration step size was necessary to achieve 
numerically stable solutions; the stability criterion is related to the hyper- 
bolic characteristic curves. 

It was found that the exact forms of the momentum and rate equations along 
the body surface could not be used in lieu of their counterparts in the approx- 
imate systems. Such a hybrid procedure always produced unstable numerical 
results near the wedge or cone tip. 


Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., June 4, 1964. 
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APPENDIX A 


FROZEN SHOCK-WAVE RELATIONS 


The appropriate frozen shock-wave relations are given for 7=7/5 
( ref . 18 ) . Let 


b(p,Mj = M,» 2 sin 2 P 


then 


= -^(b - 1) 

6Moo 


u(x,B) = Ug = (l - a) cos 0 + a cot p sin 0 


v(x,S) = Vg = -(l - a)sin 0 + a cot p cos 0 


p(x,5) 


P6 


55b - 5 


/ \ 6b 

p(x,S) = p 5 = - 

° 5 + b 

The p -derivatives of equations (A3) to (A6) are also required explicitly 
they are, respectively, 


dug 

dp 


5 „ . -v a sin 0 

- — cos p sin A - 

3 sin^p 


dv 5 5 

= — cos p cos 

dp 3 


a cos 9 
sin^p 


dp s 5 

— - = — sin p cos p 
dp 3 

^8 = 60 b cot p 
d|3 (5 + b) 2 


(Al) 

(A2) 

(A3) 

(A4) 

(A5) 

(A 6) 

and 

(A7) 

(A8) 

(A9) 

(A10) 
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APPENDIX B 


COMPUTATIONAL EQUATIONS FOR N = 2 


Before the final computational forms of the approximate systems are 
obtained, the derivatives of the functions Qi^ must be further expanded in 

terms of u-^, v^, p-^, p^, E^, and r^ . The factions at the shock wave 
Qj^g are explicit functions of p (and rg when j = 1) so that dQ-^g/dx 
can be expressed in terms of the single derivative dp/d x (and drg/dx when 
j = 1) . For example, 


d S,5 

dx 


id/ N . “o 

rS ta( P S U 5 v 6) + '^'Vs ST 


(Bla) 


or 


"£x' ~ = Sin 9 + S cos ^( p 8 u S v &)fo + j( sin 9 + tan ^ cos ®^ P 8 U 5 V 8 

(Bib) 

In any approximation the equation which contains dQl^^Qjdx and dQ^ g/dx will 

yield the required expression for dp/dx by using equations (Bib) and (ll) and 
by noting that Q^q = Pq u o v O = 0 for a11 va l ues x * Treating dp/dx ^and 

thus dQ^g/dx) as a known algebraic quantity effects considerable simplifica- 
tion^ as follows: at each strip boundary the derivatives of the corresponding 

variables u^, v^, p^-, and can be obtained by solution of a 4 by 4 linear 

algebraic system. (Eqs. ( 5 ) and (6) are used to eliminate derivatives of p k . ) 

This feature is emphasized to point out that one does not have to deal directly 
with a 4N by 4N system but rather N (4 by 4) systems. The algebra required 
for the latter is simple, so that one need not resort to time-consuming inversion 
of a 4N by 4 n matrix in the computer at each integration step. 

The procedures are now demonstrated for N = 2. The derivatives of r ^<3 
are expanded from the functions ^ ( k = 0,1,5) as. in equation (Bla), and 

then equations (l8) are divided by TqJ . Equations (l8a) can be written (in the 
order i = 3, 1,2, 4) 

SBiHp H = Kj. (B2) 
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K 




$ 

II 

'T 

0 

*1-8 

60 


(B3) 

8 to( p o + P 0 ,J 0 2 ) ■ k 5 


(B4) 

» &(e 0 V%) * K *t 


(B5) 

and equations (l8b) can be written (in the order i 

= 1,2,3,^) 


8Bs 


(b6) 

8B 2 s( p i + P 1 U 1 2 ) ■ K « 


(B7) 

8b 2 £(Wl) ‘ K 7 


(B8) 

BB2 = Kq 


(B9) 

where 



K^_ = - 2 B 2 P 1 u 1 v 1 tan A + (bj tan A - J0jpgUgVg 



- 4 [po “ B 2^Pi + Pi^i 2 ) + B^pg + PS V 5 2 ) 

- j0 cot e^p 0 - p 6 ^ 



(BIO) 

K 2 = — K]_ - (tan A + j0)p Q UQ + 2B2P 1 u 1 tan A 
H 1 

- ^Bj tan A - PgUg - 4^B 2 P 1 v 1 - BxPgVgj (Bll) 
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K 5 = p. K x - (tan A + j0)^P o + PqUq^ + 2B 2 tan * + Pi u l 2 j 

- tan A - jtf) ^> 8 + P 5 v & 2 j - 4 (b 2 P 1 u 1 v 1 - Bip 5 u 5 v & ) + j0^P o - p s j 

(B12) 

K4 = -(tan A + j0)p o u o E o + 2B 2 p 1 u 1 E 1 tan A - 4B2P 1 v 1 E 1 + S^p 0 € 0 - Bjpjgegj 

(B13) 

K5 = - ^2 Ki - ^4 tan A + 2j0^ p-jV^ + ^4 tan A - 3^jp^ + 2B 2 P 1 v 1 - 5 Bip 5 v 5 


(B14) 


Kg = - |5. K ± - ^B4 tan A + 2j0^p 1 + p-jU-^J + ^4 tan A - J0^p 5 + P 5 a 5 2 j 

+ 2B 2 p 1 u 1 v 1 - 5Bip 5 u & v s + j0^2p 1 + p 5 j (B15) 

Ky = -K]_ - ^B4 tan A + 2j0^p-^u^v^ + ^B4 tan A P & U 8 V 5 

+ p 0 + 2B 2 ^p 1 + p^^j - 5Bi^p 5 + P 5 v 5 2 ) + J0 cot e(2p ± + p 5 j (Bl6) 

Kg = -^4 tan A + 2j0)p 1 u 1 E 1 + 2B2P 1 v- L E 1 + 6^B 2 p 1 €- L + B]_p g€gj (B17) 




Hi 

h 2 

h 3 


^■( P 5 U &) 


^■(P5 + P5 U 5 2 


(B18) 


24 



B]_ = 1 + j0 cot 0^ 
B 2 = B 1 + i 

> 

B3 = 2Bi + 1 
Bl^ = 3Bi + 1 


(B19) 


The variable 0 = 6/x is introduced to maintain consistent numerical accuracy- 
in the ratio 8/x near the cone tip. 


The derivatives of 
and ( 6 ) with 7 = 7 / 5 , as 


p k u k (k = 0,1) are expanded by using equations ( 5 ) 
follows : 


dx 


(pk u k) = 




|l.4m k 2 -j— ^ + ^1 + 0.4m k 2 ^ 


du k , p 

Pk u k + °‘ 4mk l 


dv- 


* P k V k d7" + 17 


■k , p k u k 

dx 


where 


Also, 


m k = u k 



(B21) 


s( p k + vi 2 ) = Ir + p k u k !r + £( p k u k) (B22) 

^(‘Wk) - Pk“k sr + v k s( p k u k) (B23) 

£(°k"A) " p k u k ^ + % s( p k u k) ^ 


Equations (B20) to (B24) are substituted into equations (B3) to ( B 9) to arrive 
at the final equations which (together with eq. ( B2 ) ) complete the set for 
N = 2: 
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SBlHl £ = Kl 


(B2) 


6P 0 U 0 ^ = K 4 “ *0*2 


(“ o 2 " tydT " "( 1 + 1 - 4 m o 2 )( K 3 " ^ 2 ) + ” 0*2 - 5 -|- 2 - -^2 


du O d Pn k 3 - 

w sr = - dbT + -^-g — 


dEi 

5 B 2 P-JU-L — = Kg - E 1 .K 5 


dv 


U.V -I 

6 B 2 p 1 u 1 — = Ky - 


SB 2^ m l 2 ' = -(l + 1.4mi 2 )(K6 - u x K 5 j + UjKj 


P ^- v l ^1 U 1 ^1 

- °* 46B 2 m l P jY-l — - 5B 2 ^ 


du-, dp. Kg - u-tKc 

p i u i aT = " IT + 6 “ 


— = i(tan A - 0) 
dx x rJ 


(B25) 

(B26) 

(B27) 

(B28) 

(B29) 

(B30) 

(B31) 

(B32) 



APPENDIX C 


EXACT INITIAL DERIVATIVES FOR A WEDGE 


Sedney (ref, 22) used natural coordinates to derive exact expressions for 
the shock curvature and the gradients of surface flow variables at the wedge tip 
caused by vibrational relaxation. A somewhat different approach, which is more 
appropriate to the present problem, yields an equivalent solution for d(3(0)/dx 
and the initial flow-variable gradients along any strip boundary. 

The coordinates tj = y/S and £ = x are introduced to map the strips 
(fig. 2) into rectangular regions, with the shock at r\ = 1 and the surface at 
r) = 0. In the -plane the wedge tip is the line | = 0, 0 ^ r] ^ 1. Deriva- 

tives of p are eliminated from equation (l) by using equations (4), (5), and. 
(6); the transformed continuity, momentum, and vibrational-rate equations are 
as follows: 


, -vvldp , A du,dv 

uri tan A) — - t\ tan A -r— + t— 

7 P Si] St) ^ 


(JL ^ ^ + (, 

\7V d§ dg j 

s(u + (v - ut) tan A )^ 1 - T| 

( d| p dr) 


= - f (Cl) 


tan A dp _ q 
P dr) 


(C2) 


6u — + (v - ut) tan A)— + — — = 0 
d£ dr) P St) 


(C3) 


5u — + (v - ut) tan A)— = 5e 

d| dT) 


(C4) 


Equations (Cl) to (c4) are differentiated again with respect to the 
limit as | and 6 -> 0 is obtained, and the following frozen-flow wedge con- 
ditions are applied: 


v( 0 ,Tl) = 0 

on dr\ OT] 




dl) 



(C5) 


Then at the wedge tip, that is, along the line £ = 0, 
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— p 

7P\ 


" I) + ( u • "I r cot ’'I 


£ 

T 


(c6) 


pu ( u -S + p -!r 0 


(C7) 


— + pu tan a/V - t) ] = 0 
dr) \ °tW 


(c8) 


n do e 

° ^ = u 


(C9) 


where 


u - |a( 0 ,,) 
v ’ 

p = |(o,n) 

n = |f( 0 ,Ti) 


The boundary conditions are as follows: 
At T) = 1 , 


and at rj = 0. 


u(i) 

d|3 d| 

v(i) = *Sfla 

43 d| 

P(U -2££dg. 
v ' dp d| 




0 ( 1 ) = 0 


V(0) = 0 


J 


(ClOa) 


(ClOb) 
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The solution of equations (C6) to (C9) with the boundary conditions (CIO) 
is readily obtained. In terms of the notation of the main text, the wedge 
initial derivatives at x = 0 are found to be 


du k 

dx 



(j- - nQ d P& 

P 5 U 5 dp 


dp 

dx 


(Cll) 


£^k = dY S dP 

dx dp dx 


(C12) 


d Pk _ d P& dp 
dx dp dx 


(C13) 


^k 

dx 



(C14) 


dp 

dx 



$ 


1 

?8 U 6 


dpj 

W 


+ cot 


A 



-1 


(C15) 
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